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We lay out the frameworks to numerically study the structure formation in both linear and 
nonlinear regimes in general dark-matter-coupled scalar field models, and give an explicit example 
where the scalar field serves as a dynamical dark energy. Adopting parameters of the scalar field 
which yield a realistic CMB spectrum, we generate the initial conditions for our N-body simulations, 
which follow the spatial distributions of the dark matter and the scalar field by solving their equations 
of motion using the multilevel adaptive grid technique. We show that the spatial configuration of the 
scalar field tracks well the voids and clusters of dark matter. Indeed, the propagation of scalar degree 
of freedom effectively acts a fifth force on dark matter particles, whose range and magnitude are 
determined by the two model parameters (jj., 7), local dark matter density as well as the background 
value for the scalar field. The model behaves like the ACDM paradigm on scales relevant to the 
CMB spectrum, which are well beyond the probe of the local fifth force and thus not significantly 
affected by the matter-scalar coupling. On scales comparable or shorter than the range of the local 
fifth force, the fifth force is perfectly parallel to gravity and their strengths have a fixed ratio 27^ 
determined by the matter-scalar coupling, provided that the chameleon effect is weak; if on the 
other hand there is a strong chameleon effect (i.e., the scalar field almost resides at its effective 
potential minimum everywhere in the space), the fifth force indeed has suppressed effects in high 
density regions and shows no obvious correlation with gravity, which means that the dark-matter- 
scalar-field coupling is not simply equivalent to a rescaling of the gravitational constant or the mass 
of the dark matter particles. We show these spatial distributions and (lack of) correlations at typical 
redshifts {z — 0, 1, 5.5) in our multi-grid million-particle simulations. The viable parameters for the 
scalar field can be inferred on intermediate or small scales at late times from, e.g., weak lensing and 
phase space properties, while the predicted Hubble expansion and linearly simulated CMB spectrum 
are virtually indistinguishable from the standard ACDM predictions. 

PACS numbers: 04.50.Kd 



I. INTRODUCTION 

The origin and nature of dark energy Q is one of the 
most difScult challenges facing physicists and cosmolo- 
gists now. Among all the proposed models to tackle this 
problem, a scalar field is perhaps the most popular one 
up to now. The scalar field, denoted by (p, might only in- 
teract with other matter species through gravity, or have 
a coupling to normal matter and therefore producing a 
fifth force on matter particles. The latter idea has seen 
a lot of interests in recent years, in the light that such a 
coupling could alleviate the coincidence problem of dark 
energy and that it is commonly predicted by low energy 
effective theories from a fundamental theory. 

Nevertheless, if there is a coupling between the scalar 
field and baryonic particles, then stringent experimental 
constraints might be placed on the fifth force on the lat- 
ter provided that the scalar field mass is very light (which 
is needed for the dark energy). Such constraints severely 
limit the viable parameter space of the model. Different 
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ways out of the problem have been proposed, of which 
the simplest one is to have the scalar field coupling to 
dark matter only but not to standard model particles, 
therefore evading those constraints entirely. This is cer- 
tainly possible, especially because both dark matter and 
dark energy are unknown to us and they may well have 
a common origin. Another interesting possibility is to 
have the chameleon mechanism @, mil, by virtue of 
which the scalar field acquires a large mass in high density 
regions and thus the fifth force becomes undetectablly 
short-ranged, and so also evades the constraints. 

Study of the cosmological effect of a chameleon scalar 
field shows that the fifth force is so short-ranged that it 
has negligible effect in the large scale structure formation 
'6!:] for certain choices of the scalar field potential. But it 
is possible that the scalar field has a large enough mass in 
the solar system to pass any constraints, and at the same 
time has a low enough mass (thus long range forces) on 
cosmological scales, producing interesting phenomenon 
in the structure formation. This is the case of some f{R) 
gravity models 0, § , which survives solar system tests 
thanks again to the chameleon effect [l^, [ll|, ■ Note 
that the f{R) gravity model is mathematically equivalent 
to a scalar field model with matter coupling. 

No matter whether the scalar field couples with dark 
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matter only or with all matter species, it is of general in- 
terests to study its effects in cosmology, especially in the 
large scale structure formation. Indeed, at the linear per- 
turbation level there have been a lot of studies about the 
coupled scalar field and /(i?) gravity models which en- 
able us to have a much clearer picture about their behav- 
iors now. But linear perturbation studies do not conclude 
the whole story, because it is well known that the matter 
distribution at late times becomes nonlinear, making the 
behavior of the scalar field more complex and the linear 
analysis insufficient to produce accurate results to con- 
front with observations. For the latter purpose the best 
way is to perform full N-body simulations [l3| to evolve 
the individual particles step by step. 

N-body simulations for scalar field and relevant models 
have been performed before [H, [H [H, [O, [H, [M H IH ■ 
For example, in the simulation is about a specific 
coupled scalar field model. This study however does not 
obtain a full solution to the spatial configuration of the 
scalar field, but instead simplifies the simulation by as- 
suming that the scalar field's effect is to change the value 
of the gravitational constant, and presenting an justifying 
argument for such an approximation. As we will see be- 
low, this approximation is only good in certain parameter 
spaces and for certain choices of the scalar field poten- 
tial, and therefore we believe fuller simulations than the 
one performed in is needed to study the scalar field 
behavior most rigorously. Note also the structure forma- 
tion with a coupled chameleon scalar field has also been 
studied using non B-body techniques previously [22| . 

Recently there have also appeared full simulations of 
the f{R) gravity model [HjUJi which do solve the scalar 
degree of freedom explicitly. However, embedded in the 
f{R) framework there are some limitations in the gener- 
ality of these works. As a first thing, f{R) gravity model 
(no matter what the form / is) only corresponds to the 
couple scalar field models for a specific value of coupling 
strength [25I]. Second, in f{R) models the correction to 
standard general relativity is through the modification to 
the Poission equation and thus to the gravitational po- 
tential as a whole [2^, while in the coupled scalar field 
models we could clearly separate the scalar fifth force 
from gravity and analyze the former directly. Also, in 
/(i?) models the coupling between matter and the scalar 
field is universal (the same to dark matter and baryons), 
while in the couple scalar field models it is straightfor- 
ward to switch on/off the coupling to baryons and study 
the effects on baryonic and dark matter clusterings re- 
spectively. And finally, the general framework of N-body 
simulations in couple scalar field models could also han- 
dle the situation where the chameleon effect is absent and 
scalar field only couples to dark matter (which is of no 
interests for f{R) people). 

In this article we present the general formulae and algo- 
rithm of full N-body simulations in coupled scalar field 
models and consider as an explicit example the results 
for a chameleon scalar field. Unlike in [20] , here we shall 
calculate the spatial distribution of the scalar field di- 



rectly without making simplifications such as a rescaled 
gravitational constant. Neither shall we use the concept 
of varying-mass dark matter particles as in (20| , but in- 
stead we treat the system as constant-mass dark matter 
particles under the action of a fifth force. 

The article is organized as follows: in §|TT]we review the 
general equations of motion for the coupled scalar field 
model and introduce our specific choices of the coupling 
function and scalar field potential. Next we analyze in 
i^ lllll the general linear perturbation equations in the 3-1-1 
formalism [26, 27] and integrate them into the numerical 
Boltzmann code CAMB [23 to study the effects on the 
linear structure formation. Then in § IIVI we turn to our 
main focus, introducing the formulae and algorithm of 
the N-body simulation. We also study the chosen model 
explicitly, display the preliminary results and discuss on 
them. Finally we conclude in § |Vl 



II. THE COUPLED SCALAR FIELD MODEL 

In this section we briefly introduce the model consid- 
ered here and present the equations that will be analyzed 
in the following sections. Let us start by looking at the 
general field equations for a scalar field coupled to dark 
matter. 

The Lagrangian for our coupled scalar field model is 



R 

K 



V{^)-C{ip)CcDM+Csil) 



where R is the Ricci scalar, k = SttG with G Newton's 
constant, is the scalar field, V{ip) is its potential energy 
and C{ip) its coupling to dark matter, which is assumed 
to be cold and described by the Lagrangian £cdm- 
includes all the terms for photons, neutrinos and baryons, 
and these will be considered only when we calculate the 
large scale structure formation in the next section. 

The dark matter Lagrangian, for a point-like particle 
with (bare) mass mo, is 



-CDM 



(y) 



-^^S{y - xo) Jgabigij^ 



(2) 



where y is the coordinate and xo is the coordinate of the 
centre of the particle. From this equation it can be easily 
derived that 



rpah 
-'CDM 



mo 



'^(y-xo)iSio 



(3) 



Also, because gabio^^ — 9abU°^u^ = 1 where it° is the four 
velocity of the dark matter particle, so the Lagrangian 
could be rewritten as 



^CDM(y) — 

which will be used below. 



mo 



<^(y-xo), 



(4) 
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FIG. 1: Upper panel: the potential (dashed curves), coupling function (dotted curves) and effective potential (solid curves) of 
the strong coupling scalar field model at various cosmic epochs. Lower panel: the same as the upper panel but for a coupling 
not strong enough so that the scalar field will not reside at the minimum of 14//. Solid circles denote the states of the scalar 
field and in case that this does not coincide with the minimum of effective potential, the later is indicated by an open circle. 
See text for explanations. 



Eq. ([3]) is the energy momentum tensor for a single 
dark matter particle. For a fluid with many particles the 
energy momentum tensor will be 



rpab 
-'CDM 



mo 



5{y- xq)xqXq 



(5) 



in which 1^ is a volume microscopically large and macro- 
scopically small, and we have extend the 3-dimensional 
5 function to a 4-dimensional one by adding a time com- 
ponent. 

Meanwhile, using 



T 



ah 



2 H^C) 



lab 



(6) 



it is straightforward to show that the energy momentum 
tensor for the scalar field is given by 



(7) 



So the total energy momentum tensor is 



Tab = ^a^^b'-p-Qab 



T'S 



(8) 



where T^^^ = PcDM^aUb, Tf^ is the energy momen- 
tum tensor for standard model particles, and the Einstein 
equation is 



Gnh — kT„ 



(9) 



where Gab is the Einstein tensor. Note that due to the 
coupling between the scalar field if and the dark mat- 
ter, the energy momentum tensors for either will not be 
conserved, and we have 



CDMab 



CDM 



T 



CDMafcN 



Vf,(4lO) 



where throughout this paper we shall use a, to denote 
the derivative with respect to Lp. 



4 



Finally, the scalar field equation of motion (EOM) from 
the given Lagrangian is 



C. 



CDM 



dip dip 
where □ = V^Vq. Using Eq. (jH) it can be rewritten as 



(11) 



Eqs. (0 m [ini fTT|) summarize all the physics that will 
be used in our analysis. 

We will consider a special form for the scalar field po- 
tential, 



[I - exp {-y/^f)] 



(12) 



where we have fixed the coefficient in front of 93 to be 1 
without loss of generality, since we can always rescale ip 
to achieve this; is a dimensionless constant and Vb is 
a constant with mass dimension 4. As will be discussed 
below, we need <C 1 to evade local observational con- 
straints. Meanwhile, the coupling between the scalar field 
and dark matter particle is chosen as 



C{ip) = exp(7VK(^), 



where 7 > is another dimensionless constant. 

As will be explained below, the two dimensionless pa- 
rameters /I and 7 have clear physical meanings: roughly 
speaking, fj, controls the time when the effect of the scalar 
field becomes important in cosmology while 7 determines 
how important it would ultimately be. Indeed, the poten- 
tial given in Eq. (fT^ is motivated by the f{R) cosmology 
[To| . in which the extra degree of freedom behaves as a 
coupled scalar field in the Einstein frame. As we could 
see from Eq. (fT^ . the potential V ^ 00 when (p ^ 
while V Vq when — > 00. In the latter case, however, 
C ^ 00, so that the effective total potential 



^e//('/') = V{ip) + pcbmC{ip) 



(14) 



has a global minimum at some finite (p. If the total poten- 
tial Vef f (ip) is steep enough around this minimum, then 
the scalar field becomes very heavy and thus follows its 
minimum dynamically, as is in the case of the chameleon 
cosmology. If Vef / is not steep enough at the minimum, 
however, the scalar field will have a more complicated 
evolution. These two different cases can be obtained by 
choosing appropriate values of 7,/x: if 7 is very large or 
fi is small then we run into the former situation and if 7 
is small and fj, is large we have the latter. 

In Fig. [1] we present a schematic plot of the two situa- 
tions: the significant difference is that at late times when 
PCDM becomes small, the effective potential becomes flat 
around its minimum if 7 is not very large and not very 
small, and as a result the true value of ip will lag behind 
that corresponding to the minimum of Veff. Of course 



if Vq is chosen appropriately the scalar field can act as a 
dynamical dark energy in this slow-roll regime. 

The complexity of the two cases also makes them phe- 
nomenologically rich, and it is of our interests to study 
how such a setup will affect the cosmology in background, 
linear perturbation and in particular nonlinear structure 
formation regimes. In the next two sections we shall con- 
sider these questions. 



III. LINEAR STRUCTURE FORMATION 

Our first task is to linearize the model. Since this in- 
volves many equations, we refer the interested readers to 
the Appendices \X\ and [Bl for the complete first-order per- 
turbation equations and their representations in fc-space. 
We shall proceed then to study their effects on the large 
scale structure at linear perturbation level in the Universe 
by putting these linearize equations into the numerical 
Boltzmann code. 



A. The Background Evolution 



In what follows we will consider the limit of small /i, 
^^'^^ M ^ 0(0.1). It is then useful to define the effective mass 



of the scalar field by the Taylor expansion 



where ipo is at the minimum of Vef / , given by 



dVeffjcpo) 

dip 



y^ipo 



7 PCDM 



(15) 



(16) 



in which we have used the facts that 1 -I- /i « 1 as well 
as ^/Klp ^ 1 so that ex.p{—y/Kip) « eyip{'y^/Klp) « 1 and 
1 — exp{—y/K(p) w y/Kip. Thus the effective mass 



^2 



d'Veff i^o) 

difil 

nVo 
P-kVq 



y = exp{~^/Kipo) 

(i + M)y 



(17) 



1-2/ 



In- 



I'^upcDKiy ^ 



, (1 + 37)^0 



IPCDM 



IPCDM 



(18) 



To see that the scalar field is heavy, i.e., ml^j » H'^ ~ 

HQa~^, note that npcuM O (i?ga~'^) » hkVq, true for 
fi ^ 1 [or ^ 0(1) but in the early universe, remember 
that we always assume 7 ~ 0(1)]- This means that for 
/i <C 1 (or in the early times) the scalar field has a very 
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FIG. 2: (Color Online) The evolution of the fractional energy densities in radiation (red lines), baryons + cold dark matter (blue 
lines) and scalar field dark energy component (green lines) with respect to \g{l+z) where z is the redshift. The model parameters 
for this case are chosen by 7 = 1 and /i = 0.1,0.2,0.3,0.5 respectively. The physical parameters are Or = 8.475 X 10"^ for 
radiation, JIb = 0.05 for baryons, JIcdm = 0.25 for CDM and Ho — 70 km/s/Mpc for the present Hubble expansion rate. To 
ensure these parameters the value of Vo must be fine-tuned, and here we have used a trial-and-error method to adjust the value 
of A = KVo/3ii"o as 0.53285, 0.44868, 0.38968, 0.30848 respectively for fj, = 0.1, 0.2, 0.3, 0.5. Note that given Ho the propagations 
of the fractional energy densities fix the background evolution. 



heavy mass and tends to settle at ip ^ ipo (quickly oscil- 
lating there) . Then as a good approximation we have 



\ Vo J 



Vo (19) 



in which wc have used the facts that lim^^o+ /^^ = 1 
this regime. At the same time, because y^ip ~ \/Kipo ~ 
_E Vo PCDM ^ Z^m^H < H with the use of Eq. (fT6l). 

7 PCDM PCDM V -ru ^ -1 >| V^ 

so the kinetic energy of the scalar field is negligible. This 
shows that for small /i we do expect the scalar field to 
behave like a cosmological constant in background cos- 
mology. 

This analysis has been confirmed by numerical calcu- 
lations. In Fig. [5] we show the fractional energy densities 
of the radiation, dust (baryon plus CDM) and dark en- 
ergy (the scalar field) components for 7 = 1 and several 



values of It can be seen there the evolutions of these 
fractional energy densities are not sensitive to /i. To see 
this more clearly we have also plotted in Fig. [3] the total 
equation of state of all matter species. As /i ^ the 
curves quickly become indistinguishable from each other 
(and indistinguishable from the ACDM prediction). 



B. The Linear Perturbation Evolution 

To study the perturbation evolution and the effect of 
the scalar field on large scale structure formation, we just 
need to implement the perturbation equations listed in 
Appendix [B] into a Boltzmann code. As mentioned ear- 
lier, here we work in the A = frame, in which the CDM 
peculiar velocity f cdm must be dynamically evolved ac- 
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linear (on large scales) and nonlinear (on smaller scales) 
matter power spectra. 

To explain the matter power spectrum, let us consider 
the evolution of cold dark matter density contrast Acdm 
(A similar analysis is first give n in 10] for f{R) model 
and later more generally in [29|). For simplicity we shall 
assume a Universe filled with only radiation and CDM 
particles. Taking the (conformal) time derivative of the 
equation 



A' 



we have 



CDM 



A" 

^CDM 



kZ + kvi 



CDM 



kZ' + kv'c 



CDM 



= 0, 



= 0. 



(23) 



(24) 



Taking the spatial derivative of Eq. (jAlip , collecting the 
coefficients of harmonic expansions and using the frame 
choice A — we obtain 



FIG. 3: (Color Online) The total equation of state of the 
model with different values of /i as a function of lg(l + z). 
The parameters are the same as in Fig. [2] and the values of fi 
are given beneath the curves. 



kZ' 

Substituting Eq. 
we arrive at 



k-Z + -{X + 3XP)a 



, ... 0. (25) 

a 2 

5]) into Eq. ^ and using Eqs. (IB13[ 



cording to Eq. (|B13p with A = 0. Also note that the tight 
coupling approximation is not affected as compared with 
GR+ACDM. 

The initial condition for wcdm could be set to zero. 
In this way the observer is comoving with the dark mat- 
ter particles initially (when the ^ term in Eq. (|B5[) is 
extremely small) and only deviate from the dark matter 
particle world lines when this ^ term becomes significant. 
As for the initial conditions of ^ and , they are also very 
small at early times, and as we shall see below, there is 
good reason to set them to zero at those times. Thus we 
have adopted 



^cdm 



-K 



^initial 
C' 

Sinitial 
^^CDM, initial 



0, 
0, 
0, 



(20) 
(21) 
(22) 



as the initial conditions of the new variables which need 
to be evolved in the code. 

In Fig. m we show a collect of the CMB power spectra 
for the model with 7 — 0.5. In Fig. [5] the matter power 
spectra of the corresponding parameter choices are given. 
We can see that the effect of the scalar field on the CMB 
spectrum is most significant for \ow-i (large scales), and 
even there the effect is only to slightly reduce the power. 
This behavior is quite similar to the one found in 10.] and 
is due to the less decay of the gravitational potential on 
large scales at late times, which decreases the integrated 
Sachs- Wolfe (ISW) effect. Such small deviations of the 
CMB spectrum from the ACDM result are obviously not 
of much help in placing constraints on the model parame- 
ters /i and 7, especially because of the cosmic variance on 
large scales. Thus to distinguish this model from others 
we necessarily need to use other observables such as the 



CDM 



v'vcbm) = 0, (26) 



where the last term is equal to k (jtA + w^dm + ^^^cdm 
with ^ being substituted using the propagation equation 
of vcDM (remember that A = hy choice of the frame). 

Now according to Eq. (jB26P the evolution of is gov- 
erned by 

C" + 2^^' + (fc2 + a^V^^ + a^pcBuC^^)^ 

+k(p'Z + a^C^pcBM^CDM = 0. (27) 

For small scales (and late times) A; ^ ]^] and the term 
/c^^ dominates over other terms proportional to ^ in the 
above equations dominates so that the equation is ap- 
proximately 

/c^^ -I- o^C^PcdmAcdm = -kip'Z (28) 
= V^'^cdm + kip'vcDM-i29) 



where we choose to eliminate Z by using Eq. (p3|) . Sub- 
stituting Eq. ([29| into Eq. ([26|l to ehminate ^ and vcdm, 
we get 



A" 

^CDM 



I Ma' 



CDM 



c 

-|(A' + 3XP)a^ ~ -^pcdmAcdmo' = 0. (30) 

During the matter dominated era we can neglect the con- 
tribution from radiation, and the above equation for the 
growth of the overdensity reduces to 



^CDM + ^^CDM = 47rGpcDMAcDMa^, 



(31) 
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FIG. 4: (Color Online) The CMB power spectra for the couple 
scalar field model with 7 = 0.5. The red, blue, green and 
black curves correspond to ^ = 0.1,0.2,0.3 with chameleon 
perturbations and n — 0.1 without chameleon perturbations, 
respectively. Curves for /i < 10^'' are indistinguishable from 
the ACDM result. 

where 

where k, = SttG as before, and 2C^/C « 27^ for our cou- 
phng function of C as given by Eq. (I13p . Thus we could 
see that the growth of cold dark matter density contrast 
is (on small scales) scale-independent, which explains the 
matter power spectrum at large k values (small scales) 
in Fig. O In particular, note that for our choice of pa- 
rameters we have C ~ 1 and H — — + ^ ^ —, and so 

a C a ' 

this equation is approximately the same as that in the 
ACDM paradigm but with an effective gravitational con- 
stant G ~ (H- 27^ )G. This indicates that the parameter 
7 characterizes the strength of the scalar field fifth force 
relative to gravity; the larger 7 is, the stronger the fifth 
force will ultimately be. 

Meanwhile, from the derivation process of Eq. pTjl we 
see that the assumption that fc^ :s> a? [V^^ + PcdmC^^) 
[cf. Eq. [27] has been used. This is of course true for our 
choices of /i, but Eq. p7p tells us that this is not neces- 
sarily the case for ^ 1, due to the strong nonlinearity 
of V{ip). Let's suppose now fc^ < (V^^ +PcdmG^^) 
for some small enough value of /i (say /i — 10""^), then 
Eq. ([27|) should be approximated as 

{a'^V^^ + a^pcDuCipip)^ + ktp'Z + a^C^pcDW^CDM = 

rather than as Eq. 1^5]) . so that in Eq. (^5)1 the fc^^ term 
could be neglected, implying that the scalar field simply 
has negligible effects on the evolution of Aqdm- Thus 
we conclude that the effect of the parameter \i is to con- 
trol the time when the fifth force becomes important as 



FIG. 5: (Color Online) The matter power spectra for the 
couple scalar field model. The parameters are the same as in 
Fig. 131 The curves with /i < 10^^ are very close to the ACDM 
result. 

compared with gravity; the smaller /i is, the later will 
this time be. So with very strong chameleon effects (very 
nonlinear potentials) the fifth force could be greatly sup- 
pressed all through the cosmic history up to now, and it 
is therefore possible to get a cosmology very close to the 
ACDM paradigm in every aspect. 

Obviously, at the very early times the fifth force is also 
suppressed even in cases of p ^ C(l), because in this 
regime again we have <C (VJ^i^ -I- pcDiviGipi^). Indeed 
this is the reason why we could set the initial condition 
as in Eq. pl)l . 

Now the matter power spectra as seen in Fig. [5] could 
be well explained: when p decreases, the time when the 
fifth force could fully show its power becomes later and as 
a result the growth of structure does not reach its max- 
imum potential - this explains why the /i = 0.1,0.2,0.3 
curves increasingly get far away from the black one. But 
because all these values of p are large enough (compared 
with, say, \i = 10""^), the fifth force in all these cases has 
almost realized its full power (which is 27^ times stronger 
than gravity according to Eq. pip ) - this explains why 
all the colored curves are fairly far from the black one. 
Finally, Eq. ((3T|) explains why on small scales the colored 
P(fc) curves are almost parallel to the black one. 

Given the complexity of the scalar field behavior, es- 
pecially in regions with highly nonlinear matter distri- 
bution, we must say that the linear analysis performed 
above is only qualitatively correct and lacks the high pre- 
cision to compare with various cosmological data sets. In 
particular, it may well be that the effect of the fifth force 
is not significant on linear scales because the fifth force 
is not long-range enough, but much more important on 
smaller scales which are beyond the linear regime. Thus 
in the next section we shall consider the scalar field model 
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in the context of N-body cosmological simulations, and 
try to study its effects in a more accurate manner. 



IV. NONLINEAR STRUCTURE FORMATION 

In § mil we have considered the background evolution 
and linear large structure formation in the couple scalar 
field model. For the former, the coupling between Lp and 
CDM prevents (f from rolling to infinity and instead tends 
to keep Lp and V{}p) constant; here we find that for a value 
of /i < 0(0. 1) this effect is significant enough to make the 
background cosmology similar to that of ACDM. For the 
latter, the CMB spectrum is also not very sensitive to 
the value of fi or 7 such that < 0(0.1) is difficult to 
be distinguished from ACDM. These results suggest that 
large scale observables are not particularly well suited in 
studying the new features of this class of coupled scalar 
field models. 

On smaller but still linear scales, we already see that 
the scalar field coupling effectively increases the gravita- 
tional force by a factor of 27^, leading to a significant in- 
crease in the small scale power of P{k) for gravitational- 
strength couplings [7 ~ 0(1)]. One could of course argue 
that the observed matter power spectrum is indeed for 
the luminous matter only and cannot be applied to dark 
matter naively due to the bias between the two. But the 
enhancement of dark matter clustering in this model is 
abstract and can be tested by observations such as weak 
lensing without much ambiguity. If this is done, it is un- 
likely that there remains much space for such a strongly 
coupled model. 

However, two things must be taken into account before 
we arrive at any definite conclusions about the fate of the 
model. First, on small scales where the scalar field effect 
becomes important, the distribution of matter is already 
beyond the linear perturbation regime and in some re- 
gions could be very nonlinear. This means that a linear 
analysis as presented in the above section is no longer suf- 
ficient and we must consider the nonlinear effect, which is 
most precisely taken in account by large N-body simula- 
tions. This will be the very topic of this section. Second, 
as we argued in § IIIIl the behavior of the model is not 
only controlled by 7, but also by parameter /x: for /i <C 1, 
the epoch when the scalar field fifth force starts to real- 
izes its full power as 27^ times of gravity could be greatly 
postponed. It is then possible to prevent too powerful a 
structure formation. 

Like in the background cosmology, this has something 
to do with nonlinear (chameleon) effects, but here things 
become much more complicated. In a homogeneous uni- 
verse, the spatial gradient of the scalar field vanishes and 
the scalar field takes the same value anywhere, largely 
simplifying the analysis. For the real universe with non- 
linear matter distributions, the spatial gradient of (p is 
normally much larger than the time derivative so that 
the configuration of the scalar field relies on the under- 
lying matter distribution sensitively in a nonlinear way. 



and at the same time strongly affects the latter via the 
action of fifth force. Quantifying such complex couplings 
also calls for the use of N-body simulations. 

In this paper we shall set up the basic framework of 
studying coupled scalar field models using N-body tech- 
nique, putting much emphasis on the working mechanism 
and the fifth force effect. We shall start with a derivation 
of all the relevant equations of motion, describe how to 
integrate them into the numerical code and present some 
preliminary results. Detailed analysis will be postponed 
to companion papers. We will not provide an introduc- 
tion to general N-body simulation techniques here, and 
interested readers are referred to the relevant literature 



A. The Nonrelativistic Equations 

Our first step is to simplify the relevant equations of 
motion in the appropriate limit to get a set of equations 
which can be directly applied to the numerical code. 

We know that the existence of the scalar field and its 
coupling to standard cold dark matter particles make the 
following changes to the ACDM model: First, the energy 
momentum tensor has a new piece of contribution from 
the scalar field; second, the energy density of dark mat- 
ter is multiplied by a function C((^), which is because the 
coupling to scalar field essentially renormalizes the mass 
of dark matter particles; third and most important, dark 
matter particles will not follow geodesies in their motions 
as in ACDM, rather, the total force on them has a con- 
tribution from the scalar field. 

These imply that the following things need to be mod- 
ified or added: 

1. The scalar field Lp equation of motion, which deter- 
mines the value of the scalar field at any given time 
and position. 

2. The Poisson equation, which determines the grav- 
itational potential (and thus gravity) at any given 
time and position, according to the local energy 
density and pressure, which include the contribu- 
tion from the scalar field (as obtained from ip equa- 
tion of motion) now. 

3. The total force on the dark matter particles, which 
is determined by the spatial configuration of ip, just 
like gravity is determined by the spatial configura- 
tion of the gravitational potential. 

We shall describe these one by one now. 

For the scalar field equation of motion, we denote tp as 
the background value of Lp and 5lp = Lp — (p as the scalar 
field perturbation. Then Eq. (Ilip could be rewritten as 

5lp + ZHSlp + vIlp + V^{lp) - V^{lp) 
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by subtracting the corresponding background equation 
from it. Here V^a is the covariant spatial derivative with 
respect to the physical coordinate r = ax with x the con- 
formal coordinate, and = VraV"- Vra is essentially 
the Va, but because here we are working in the weak field 
limit we approximate it as = — ^9^^ + d"^^ + d^^ ^ by 
assuming a flat background; the minus sign is because our 
metric convention is (+, — , — , — ) instead of (— , +, +, +). 
For the simulation here we will also work in the quasi- 
static limit, assuming that the spatial gradient is much 
larger than the time derivative, |Vr(/3| ^1^1 (which will 
be justified below). Thus the above equation can be fur- 
ther simplified as 



(33) 



in which 9^ = = + {dx + dy+ d^) is with respect 

to the conformal coordinate x so that Vx = aVr, and 
we have restored the factor in front of (the Lp here 
and in the remaining of this paper is times the tp in 
the original Lagrangian unless otherwise stated). Note 
that here V and pcDM both have the dimension oimass 
density rather than energy density. 

Next look at the Poisson equation, which is obtained 
from the Einstein equation in weak-field and slow-motion 
limits. Here the metric could be written as 

ds^ = (1 + '2.(j))dt^ - (1 - 2ijj)5,jdr^dr^ (34) 

from which wc find that the time-time component of the 
Ricci curvature tensor = — V^0, and then the Ein- 
stein equation Rab = n {Tab ~ \gabT) gives 

i?°o = -Vr0 = -^(ptot + 3ptot) (35) 

where ptot and ptot are respectively the total energy 
density and pressure. The quantity Vrf/) can be expressed 
in terms of the comoving coordinate x as 



^16 = Ivi \aa^^ 



1 



a 2 



a 



= — Vi$-3- 
where we have defined a new Newtonian potential 

^ = a(j)+ ifl^ax^ 
and used Vjx^ = 6. Thus 



(36) 



(37) 



where in the second step we have used Eq. ([35)1 and the 
Raychaudhrui equation, and an overbar means the back- 
ground value of a quantity. Because the energy momen- 
tum tensor for the scalar field is given by Eq. ([7]), it is 
easy to show that p*^ -I- Sp"^ — 2 [i^^ — V^((p)] and so 

= -4TrGa^pcDMCip>) + 2[ip^ -V{p)]} 
+47rGa3 {pcDMC((p) + 2[f- Vi^)] } . 



Now in this equation (p^ — (fP' = 2(p6ip + Sip <C (Vrfp)^ 
in the quasi-static limit and so could be dropped safely. 
So we finally have 



= AirGa^ [pcBMC{ip) - pcBMC{ip)] 
-SttGo^ [V{ip) - V{ip>)] . 



(39) 



Finally, for the equation of motion of the dark matter 
particle, consider Eq. pO|) . Using Eqs. ([3l[4]), this can be 
reduced to 
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b\ CVMt7 ,n 



(40) 



Obviously the left hand side is the conventional geodesic 
equation and the right hand side is the new fifth force due 
to the coupling to the scalar field. Before going on fur- 



ther, note that the fifth force 



Civ) 



Vaf = ValogC(^) 



is perpendicular to the 4- velocity Ua] this means that 
the energy density of CDM will be conserved and only 
the particle trajectories are modified as mentioned above. 
Now in the non-relativistic limit the spatial components 
of Eq. (|40p can be written as 



(41) 



where t is the physical time coordinate. If we instead use 
the comoving coordinate x, then this becomes 



2-x = -^Vx* 



1 c^iv) 

a2 dip) 



Vx(p (42) 



where we have used Eq. (j37p . The canonical momentum 
conjugate to x is p = a^x so we have now 



dx 
dp 



a 



C{ip) 



V, 



(43) 
(44) 



Eqs. (131 [Ml m HI) wiU be used in the code to evaluate 
the forces on the dark matter particles and evolve their 
positions and momenta in time. 



(38) 



2 (ptot + 3ptot) — 2 (ptot + 3ptot) 



B. Internal Units 

In our numerical simulation we use a modified version 
of MLAPM (^3fli], see IIV Cp . and we must change our 
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above equations in accordance with the internal units 
used in that code. Here we briefly summarize the main 
features. 

MLAPM code uses the following internal units (with 
subscript c)- 

Xc = x/B, 
p, = p/{HoB) 
tc — tHo 

Pc = p/p, (45) 

in which B is the present size of the simulation box and 
Ho is the present Hubble constant. Using these newly- 
defined quantities, it is easy to check that Eqs. 
[39l [33l) could be rewritten as 



dxc _ Pc 
dtc ' 

ate o, C 



C 



(46) 
(47) 
^a', (48) 



and 



{BHo) 
3 



2 (a<p) 



C 



1 



^^a^ (49) 



where ficDM is the present CDM fractional energy den- 
sity, we have again restored the factor and again the (f 
is times the (p in the original Lagrangian. Also note 
that from here on we shall use V = d^^, = d^^ ■ d^^ 
unless otherwise stated, for simplicity. 
We also define 



X 
u 

Vo 



ln(e^ - 



1) 



(50) 
(51) 
(52) 



to be used below. 

Making discrctized version of the above equations for 
N-body simulations is non-trivial task. For example, the 
use of variable u instead of ip (see below) helps to pre- 
vent iy9 < 0, which is unphysical, but numerically possible 
due to discretization. We refer the interested readers to 
Appendix [Cl to the whole treatment, with which we can 
now proceed to do N-body runs. 



C. The N-Body Code 

The full name of MLAPM is Multi-Level Adaptive 
Particle Mesh code. As the name has suggested, this 



code uses multilevel grids [3l|, |33, l33| to accelerate the 
convergence of the (nonlinear) Gauss-Seidel relaxation 
method [11] in solving boundary value partial difi^eren- 
tial equations. But more than this, the code is also 
adaptive, always refining the grid in regions where the 
mass/particle density exceeds a certain threshold. Each 
refinement level form a finer grid which the particles will 
be then (re)linked onto and where the field equations will 
be solved (with a smaller time step). Thus MLAPM has 
two kinds of grids: the domain grid which is fixed at 
the beginning of a simulation, and refined grids which 
are generated according to the particle distribution and 
which are destroyed after a complete time step. 

One benefit of such a setup is that in low density re- 
gions where the resolution requirement is not high, less 
time steps are needed, while the majority of computing 
sources could be used in those few high density regions 
where high resolution is needed to ensure precision. 

Some technical issues must be taken care of however. 
For example, once a refined grid is created, the particles 
in that region will be linked onto it and densities on it 
are calculated, then the coarse-grid values of the grav- 
itational potential are interpolated to obtain the corre- 
sponding values on the finer grid. When the Gauss-Seidel 
iteration is performed on refined grids, the gravitational 
potential on the boundary nodes are kept constant and 
only those on the interior nodes are updated according to 
Eq. (jClip : just to ensure consistency between coarse and 
refined grids. This point is also important in the scalar 
field simulation because, like the gravitational potential, 
the scalar field value is also evaluated on and communi- 
cated between multi-grids (note in particular that differ- 
ent boundary conditions lead to different solutions to the 
scalar field equation of motion). 

In our simulation the domain grid (the finest grid that 
is not a refined grid) has 128'^ nodes, and there are a 
ladder of coarser grids with 64^, 32^, 16^, 8^, 4'^ nodes 
respectively. These grids are used for the multi-grid ac- 
celeration of convergence: for the Gauss-Seidel relaxation 
method, the convergence rate is high upon the first sev- 
eral iterations, but quickly becomes very slow then; this 
is because the convergence is only efficient for the high 
frequency (short-range) Fourier modes, while for low fre- 
quency (long-range) modes more iterations just do not 
help much. To accelerate the solution process, one then 
switches to the next coarser grid for which the low fre- 
quency modes of the finer grid are actually high frequency 
ones and thus converge fast. The MLAPM solver adopts 
the self-adaptive scheme: if convergence is achieved on a 
grid, then interpolate the relevant quantities back to the 
finer grid (provided that the latter is not on the refine- 
ments) and solve the equation there again; if convergence 
becomes slow on a grid, then go to the next coarser grid. 
This way it goes indefinitely except when converged solu- 
tion on the domain grid is obtained or when one arrives at 
the coarsest grid (normally with 2^ nodes) on which the 
equations can be solved exactly using other techniques. 
For our scalar field model, the equations are difficult to 
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solve anyway, and so we truncate the coarser-grid series 
at the 4^-node one, on which we simply iterate until con- 
vergence is achieved. 

For the refined grids the method is different: here one 
just iterate Eq. (jCll|) until convergence, without resort- 
ing to coarser grids for acceleration. 

As is normal in the Gauss-Seidel relaxation method, 
convergence is deemed to be achieved when the numerical 
solution after n iterations on grid k satisfies that the 
norm || • || (mean or maximum value on a grid) of the 
residual 

e'= = lHu';,) - h, (53) 
is smaller than the norm of the truncation error 

r^- = L'-\nu':,)^n[L\u^^)] (54) 

by a certain amount. Note here L'^ is the discretization 
of the differential operator Eq. (|C9|) on grid k and L*^^^ 
a similar discretization on grid k — I, fk is the source 
term, TZ is the restriction operator to interpolate values 
from the grid k to the grid k — 1. In the modified code 
we have used the full-weighting restriction for TZ. Corre- 
spondingly there is a prolongation operator V to obtain 
values from grid k — 1 to grid k, and we use a bilinear 
interpolation for it. For more details see [soj . 

MLAPM calculates the gravitational forces on parti- 
cles by centered difference of the potential $ and propa- 
gate the forces to locations of particles by the so-called 
triangular-shaped-cloud (TSC) scheme to ensure momen- 
tum conservation on all grids. The TSC scheme is also 
used in the density assignment given the particle distri- 
bution. 

The main modifications to the MLAPM code for our 
model are: 

1. We have added a parallel solver for the scalar field 
based on Eq. (jC3p . The solver uses a similar non- 
linear Gauss-Seidel method and same criterion for 
convergence as the Poisson solver. 

2. The solved value of u is then used to calculate lo- 
cal mass density and thus the source term for the 
Poisson equation, which is solved using fast Fourier 
transform. 

3. The fifth force is obtained by differentiating the u 
just like the calculation of gravity. 

4. The momenta and positions of particles are then 
updated taking in account of both gravity and the 
fifth force. 

There are a lot of additions and modifications to ensure 
smooth interface and the newly added data structures. 
For the output, as there are multilevel grids all of which 
host particles, the composite grid is inhomogeneous and 
thus we choose to output the positions, momenta of the 
particles, plus the gravity, fifth force and scalar field value 
at the positions of these particles. We can of course easily 
read these data into the code, calculate the corresponding 
quantities on each grid and output them if needed. 



D. Preliminary Numerical Results 

In this subsection we shall present some preliminary 
results of several runs and give a sense about the quali- 
tative behaviors of the coupled scalar field model. Also 
we will not make detailed and quantitative analysis in 
this paper, which will be shown in forthcoming papers. 

We have performed 8 runs of the modified code with 
parameters 7 = 0.5,1 and = IQ-"^, 10-^ 10^6, 10"^ 
respectively. For all these runs there are 128'^ dark matter 
particles, the simulation box has a size B = 64:h~^ Mpc 
in which h = i/o/(100 km/s/Mpc) and 128 domain grid 
cells in each direction. We assume a ACDM background 
cosmology which is a very good approximation for fi <ti 1 
as we mentioned in § IIIII the current fractional energy 
densities of dark matter and dark energy are Hcdm = 
0.28 and Ha — 0.72 (note that this is a dark-matter-only 
simulation and baryons will be added in a later work to 
study the bias effect caused by the dark matter coupling). 

Given these parameters, the mass resolution of the sim- 
ulation is 9.71 X 10^ Mq with Mq the solar mass. The 
spatial resolution is ^ 23AAh~^ kpc on the finest refined 
grids and 0.5h~^ Mpc on the domain grid. The high res- 
olution in high density regions is actually necessary to 
ensure precision in those regions because the fifth force 
is generally short-ranged there. 

All the simulations start at redshift 2; = 49. In princi- 
ple, modified initial conditions (initial displacements and 
velocities of particles which is obtained given a linear 
matter power spectrum) need to be generated for the cou- 
pled scalar field model, because the Zel'dovich approxi- 
mation [33, [ill is also affected by the scalar field coupling. 
In practice, however, we find that the effect on the lin- 
ear matter power spectrum is negligible (< C(10~^)) for 
our choices of parameters 7, /i. Thus we simply use the 
ACDM initial displacements/velocities for the particles 
in these simulations, which are generated using GRAFIC 
[3d] again with I^cdm = 0.28 and Ha = 0.72. as = 0.88 
at present day. 

In Figs, m [7] we have shown the results for the runs 
with 7 = 0.5,1.0 and /i — 10^^,/j,~^. These choices of 
/i are such that for n = 10~^ (the lower two rows) the 
chameleon effect is pretty strong while for /i = 10^^ (the 
upper two rows) it is much weaker; the choices of 7 are 
to see the effects of different full strengths of the fifth 
force. The three panels in the first and third rows display 
the particle distributions at three output redshifts z = 
5.5, 1.0, 0.0 from right to left, and the three panels in 
the second and fourth rows show the correlation between 
the fifth force and gravity at these redshifts. 

For clearness we have only plotted a thin slice (along 
the X direction) of the full 3-dimensional particle distri- 
bution. First let's have a look at the first and third rows. 
Because we have output the scalar field value at the po- 
sition of each particle together with other parameters of 
the particle, we also include this information in the plots. 
On top of each panel the range of the value C{(p) at the 
positions of all the shown particles is shown, and color 
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FIG. 6: (Color Online) The simulation results for the 7 = 0.5 and fj, = 10~^, 10~^ runs at redshifts z — 0, z = 1, and z = 5.5. 
Shown are the particle spatial yz distribution in a slab of a; = 31.5 — 32.5Mpc/h (upper panel) and the Ig-lg diagram of the fifth 
force vs. gravity in this slab (lower panel). See text for a detailed description. Each particle is a symbol with its color denoting 
the the value of the effective mass C(ip), whose minimum/maximimu are given on top of each panel, and correspond to the 
two ends of the color scale shown on the top. The size of square symbols in lower panels are proportional to the mis-alignment 
angle between the fifth force and the gravity on a particle; the biggest squares correspond to anti-alignment, and particles 
with well-aligned forces are shown as dots. A line showing a fifth-to-gravity ratio of 27^ is also drawn to show the (lack of) 
correlations of the two forces. Note that C{ip) is generally the biggest in voids where the forces are the weakest, poorly aligned 
and less-correlated. 



13 



iU-.OOO 1<G< 1.0CH>02a7 

SO 

aa 

10'' 

L...- J ■!i^...^.U-:^.J:l— Li--J.fci 



D 10 So JO *CI iWHiKA 

gainmi3=l-0ml*-5 




iO.DOD li<C< l.OOODOai' 




i 1.029 l<d< 1. 000004 7 




D 10, So JO *0 KMpcA' 

^mfflo= 1 .Om l*-5 




.4 -J -2-19 1 



25.539 t<C< 1. 0000002 

6Q 
50 
■fO 
10 

20 

10 




10 in. lit 40 MWpe/h 
ga[yifTii3= 1 -Om l*-5 




-fl.D -3,i-J.O-2,5 -2,0 -1J 
Lg(9nr) 



i 1,010 1<C< 1.00O0001 



60 
bO 

JO 

20 



laie 





10 So JO *0 iO-fjx/h 



10, So *a sfWiKAi 



25.539 t<C< I.OOOOOOO 

60 ■ " ■ " ■ ■■ 

50 
.40 
JO 
20 

10 

IL. :.L:^..jLi!T,..it^ 



10 so, J6 40 MMpc/h 




Lg{grv) Lg(grvJ Lg(grv) 



FIG. 7: (Color Online) The same as Fig.[6]but for the 7 = 1.0 and fj, = 10"^ 10"^ runs. 
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is used to illustrate the amplitude of C{lp) (going from 
black to red from the minimum to the maximum value 
of C{(p)). Note that C{lp) — exp{'jy/K(p) = 1 + jy^(p 
for ^/k(p ^ 1, thus the color also indicates the value of 
the scalar field indirectly. Also bear in mind that the 
same color may denote different values of C{lp) at differ- 
ent redshifts. The scalar field or C{(p) grows with time 
in general. 

Next look at the second and fourth rows, which display 
the logarithmic of the magnitude of the fifth force versus 
that of gravity. The color here has the same meaning as 
above. Furthermore, each point (particle) now is engaged 
with a square box centered on it, which denotes the size 
of the angle (from to tt) between the two forces; the 
size of the box increases linearly with the angle, from a 
minimum size to a maximum size comparable to the 
largest box size shown in the entire figure. As we have 
mentioned above, the strength of the fifth force, if not 
suppressed by the chameleon mechanism, is 27^ times 
that of gravity, so we also plot the functions 

IgF - lgG + lg(272), (55) 

where F, G are respectively the magnitudes of the fifth 
force and gravity, as the black solid lines in these figure, 
to compare with the simulation data. 

We can understand these results qualitatively as fol- 
lows by taking Fig.[S]as example. The chameleon effect is 
generally stronger at earlier times when matter density is 
high and the scalar field value is small. As is shown in the 
panels of the first row, at redshift z — 5.5 there is a strong 
contrast of the value C'(ip) in high density (blue) regions 
(clusters hereafter) compared to C{ip) in the low density 
{green, yellow and red) regions (voids hereafter), which 
is a direct reflection of the nonlinearity in the scalar field. 
As time evolves and the background value of the scalar 
field increases, at redshift z = 1.0 the chameleon effect 
gets suppressed; this is manifested by the facts that (1) in 
the regions of small clusters the scalar field values are no 
longer significantly different from the background value, 
both of which are yellow and orange-colored, (2) even in 
the largest clusters the contrast between the scalar field 
values inside and outside (green/light blue versus yel- 
low/orange) is not so strong compared with the result at 
redshift z = 5.5 (purple/dark blue versus green/yellow). 
These show that the size of nonlinear regions is shrinking 
and the thin shells in the clusters are thickening, together 
leading to less nonlinear behaviors of the scalar field. At 
redshift z = 0.0 this tendency just becomes more obvious, 
leaving only a small portion of the space with chameleon 
effect and thin shells. 

The strength of the chameleon effect is determined 
by several factors. In principle, the larger the effective 
mass of scalar field (mg//) at a position is, the shorter- 
ranged the fifth force is ^ and the less sensitive the scalar 
field value at this position will be to the matter distri- 
bution around it: this in turn corresponds to a stronger 
chameleon effect because the scalar field value is mainly 
determined by the local matter density. On the other 



hand, me// depends on (smaller fj, implies larger TOg// 
within a given cluster and thus stronger chameleon ef- 
fect), 7 and local pcDM (larger values for these two pa- 
rameters also imply larger m^f / and stronger chameleon 
effect) and also the background value of the scalar field 
(this is because the interior solution of the scalar field 
inside a cluster should only be solved given the bound- 
ary conditions outside the cluster, due to the differential 
nature of the scalar field EOM. Of course for very small 
fi and/or very large 7, the scalar field EOM indeed be- 
haves as an algebraic equation and then the influence of 
the background value of the scalar field is not important, 
but for our choices of ^,7 that influence is significant). 
These analysis agree well with what we have seen in the 
scalar field configuration from the above N-body simula- 
tion results (first and third rows of Fig.[S|). Also note that 
the fifth force is much weaker in regions where chameleon 
effect is strong, because a particle there can only feel the 
fifth force from those particles that are very close to it: in 
this sense we say that the chameleon effect could suppress 
the fifth force. 

A comparison between gravity and the fifth force can il- 
lustrate this more clearly. Remember that the chameleon 
effect could strongly suppress the fifth force. At low red- 
shift {z — 0.0) the chameleon effect is not significant and 
so the fifth force is not suppressed; in this case we find 
that there is a strong correlation between both the mag- 
nitudes and directions of these two forces, and the simu- 
lation results agree with the prediction Eq. (|55p to very 
high degrees. Both forces are stronger in clusters and 
weaker in voids as expected. Going backwards in time to 
the redshift z — 1.0, there is still a good correlation but 
the simulation results begin to scatter over the predicted 
line Eq. in the void regions due to the chameleon ef- 
fect (the scalar field potential term in Eq. becomes 
comparable with the matter coupling term). Finally, at 
very early times (z = 5.5) the scalar field value becomes 
so small that the potential term in Eq. (1^5]) is indeed 
much more important than the matter coupling term so 
that the latter can be neglected: we then have a complete 
mismatch between the numerical results and the predic- 
tion Eq. ([55]) as almost all data points are significantly 
below the black solid line. 

The above mismatch in the strongly chameleon regime 
can be understood schematically as follows: 

1. When there is no potential for the scalar field but 
just a matter coupling, the fifth force is indeed long 
ranged and can probe the same region as gravity 
does. Then a comparison of Eqs. (|T7l 115111^ shows 



^ Note that strictly speaking the fifth force between two particles in 
this model depends on the detailed matter distribution between 
these particles, and it is not very accurate to simply relate it to 
a certain scalar field mass m^ff. However, the concept of a fifth 
force whose range is ~ ^eff qualitatively correct and can help 
understand the situation more intuitively. 
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that the fifth force is exactly 27^ times of gravity. 
This is indeed what we have observed for z = 0.0 
when the potential term in Eq. (j49|) is negligible. 

2. When the scalar field has a potential, it acquires an 
effective mass which is well-known to be inversely 
proportional to the range of the scalar fifth force. 
As a result the fifth force is no longer as long range 
as gravity. Meanwhile, Eq. (jTZ]) makes it clear that 
the scalar field (or equivalently ln[C((y5)]) acts as 
a potential for the fifth force, and Eq. (|49p shows 
that this potential depends on the underlying mat- 
ter distribution in a different (nonlinear) way from 
what the gravitational potential <i> does. 

So we could see that there will be differences be- 
tween both the ranges and the magnitudes of grav- 
ity and fifth force. In high density regions these dif- 
ferences will be dominated over by the competing 
effects that (i) the majority part of the (either grav- 
itational or fifth) force on a particle is contributed 
by nearby particles and there are so many parti- 
cles nearby that contribution from distant parti- 
cles is negligible, (ii) in Eq. (|49p the potential term 
is much less important than the matter coupling 
term, together making the fifth force behavior sim- 
ilar to that of gravity again. In the void regions the 
number of nearby particles is small so that contri- 
bution from distant particles should be taken into 
account, and the matter coupling term in Eq. (|49|) 
is much smaller, then the difference between grav- 
ity and the fifth force becomes manifesting. These 
effects can be observed in the z = 1.0 panel. 

3. At even higher redshifts the potential of the scalar 
field makes its effective mass very large and thus 
its range very short compared with that of gravity. 
Meanwhile Eq. (|49)) becomes very nonlinear due to 
the smallness of the scalar field value, and thus the 
fifth force potential depends on matter distribution 
very differently from the gravitational potential. 

As a result, the fifth force reflects the matter distri- 
bution in a very small region around a given parti- 
cle, while gravity probes that in a much larger re- 
gion. Even in that small region where both forces 
exist, their magnitude can generally be very differ- 
ent. So it is not surprising that the two forces look 
so different as in the z = 5.5 panel. Note that here 
the fifth force is also much weaker than gravity be- 
cause its strength is suppressed by the smallness of 
the scalar field value. 

Result for the /i = 10^^ case (the lower two rows) is 
qualitatively the same as that of the = 10"^ case above, 
but here because ^ is much smaller, so the chameleon ef- 
fect exists until much more recent than in the ^ = 10""" 
case. Indeed, even at low redshifts z = 1.0 and 0.0 there 
is still strong contrast between the scalar field values in- 
side and outside the clusters (purple/dark blue versus 
yellow/orange). The correlation between the forces also 



follow our above analysis, but here there are some new 
features. The first feature is that even today the fitting 
to Eq. ([55]) is far from perfect; this is easy to understand, 
because the scalar field potential is so nonlinear that the 
scalar field potential term in Eq. (|49p is important up 
to now. The second feature is that in the z — 0.0 (also 
z = 1.0) panel we could find that in high density regions 
the fifth force does not obey Eq. (|55|) as in the fi = 10^^ 
case, but becomes much smaller than gravity - actually, 
the stronger gravity is, the weaker the fifth force will be! 
This is again due to the strong chameleon effect in the 
clusters, which makes the scalar field value very small 
and thus suppresses the fifth force there. Note that this 
feature is desirable because if we also couple the scalar 
field to baryonic matter then we definitely want the fifth 
force to be suppressed to evade solar system tests. 

If we choose 7 = 1.0 as in Fig. [3 then all the quali- 
tative results we have obtained above should still apply. 
But the stronger coupling between matter and the scalar 
field enhances the chameleon effect. This implies that a 
coupling whose strength is significantly larger than that 
of gravity could probably produce the correct amount of 
large scale structure as observed by suppressing the fifth 
force on all cosmological epochs and scales of interests to 
us. 

The analysis above clearly shows the complexity of the 
fifth force and its possible effects on the nonlinear struc- 
ture formation. Though in certain (no chameleon) limits 
the fifth force is greatly simplified and one can assume 
a modified gravitational constant in the N-body simula- 
tions as an approximation, this is evidently not the case 
if the scalar field potential is too much nonlinear. Of 
course, we have no a priori knowledge about when the 
chameleon effect becomes important, and full numerical 
simulations like the one presented here are therefore nec- 
essary to precisely study the effects of (coupled) scalar 
fields in the large scale structure. In forthcoming works 
we shall analyze the nonlinear matter power spectrum, 
halo profile, scalar field configuration within clusters, as 
well as the development and disappearance of thin shells 
in a quantitative manner. 

One may wonder if a more complete simulation should 
keep the time derivatives of the scalar field in Eq. ([55]) . 
This is certainly true, yet these terms indeed have negligi- 
ble effect [231 , which could be understood in the following 
way: when there is no (or very weak) chameleon effect, 
the scalar field potential term in Eq. is negligible 
and so this equation has the same form as the Poisson 
equation Eq. (j48p : as a result the quasi-static approxi- 
mation works as well for the scalar field as for the gravi- 
tational potential (which is what any N-body code relies 
on). On the other hand, if there is strong chameleon ef- 
fect, the scalar field value tends to be much smaller which 
means that the time derivative of the scalar field also gets 
smaller; and at the same time the spatial gradient of the 
scalar field becomes larger: these together indicate that 
the quasi-static approximation should be good here too. 
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V. CONCLUSION 

To conclude, in this paper we have presented the gen- 
eral frameworks to study the linear and nonlinear struc- 
ture formations in coupled scalar field models, and given 
some preliminary numerical results for both in the con- 
text of a specific coupling function and scalar field po- 
tential. 

For the linear large scale structure, we write down the 
perturbed field equations using the 3-1-1 decomposition, 
which can be directly applied into numerical Boltzmann 
codes such as CAMB to generate the CMB and matter 
power spectra. For the chosen coupling function with 
parameter 7 and potential with parameter /Lt, we find 
that 7 roughly controls the strength of the fifth force 
(which is due to the propagation of the scalar field) while 
/Lt controls how much the effect of the fifth force is sup- 
pressed by the chameleon mechanism. With fj, < O(O.l) 
the chameleon effect makes the model behave like ACDM 
on very large scales, but this is far from enough to signif- 
icantly decrease the effects of the fifth force on the small 
scale density perturbation growth. On those small scales, 
however, nonlinearity becomes an important issue, which 
leads us to the N-body simulation in § IIVI 

Previous N-body simulations with scalar fields are gen- 
erally simplified by certain approximations such as treat- 
ing the scalar coupling effect as a simple change of grav- 
itational constant G, or assuming a Yukawa-type force 
with a certain range. These approximations do not work 
well in the present model, and so here we set up the 
formula needed for a more precise simulation, putting 
much emphasis on the calculation of the fifth force and 
its action on particles. The equations derived and the al- 
gorithm described in § IIVI are general enough and should 
be directly applicable to other coupled scalar field mod- 
els. 

We integrate the equations into a modified version of 
the N-body code MLAPM and performed several runs 
with different combinations of 7, /j,. Some results are dis- 
played in Figs. [H] and [T] It is confirmed that when the 
chameleon effect is not important, the fifth force is par- 
allel to gravity and is 27^ times stronger, meaning that 
simply using a different gravitational constant which is 
27^ -I- 1 times as large as the bare one in the simulation 
should be an acceptable approximation. There are some 
caveats however. For one thing, the correlation between 
gravity and the fifth force is only good for a portion of 



the parameter space (7,/x), and for // <C 1 the poten- 
tial is so nonlinear as to destroy this correlation. More 
importantly, even the correlation is perfect now it is pos- 
sible that at earlier times the nonlinear effects modify it 
dramatically (cf. Fig. [7] with ^ = 10"^ at z = 0.0,1.0). 
This means that the above-mentioned approximation is 
unlikely to be correct consistently and full simulations as 
the ones in this paper are needed. 

All in all, from the results we can spot the trend that, 
with the same value of fi increasing 7 simply enhances the 
power of structure growth, while with the same value of 
7 decreasing fi has the effects of reducing that power. 
Also, for small /x the configuration of the scalar field 
becomes very nonlinear and highly sensitive to the un- 
derlying matter distribution, while for large fi this be- 
comes much more smooth and stiff. All these observa- 
tions (and others as explained in § lIVPp agree with the 
chameleon analysis. Our simulations also provides a test 
bed for scalar field theories. Among the applications, one 
could look for high-speed encounters of halos to see the 

fi^robability of generating a Bullet-cluster-like encounters 
16, 37]. Since this paper is only served to set up the gen- 
eral framework of linear and nonlinear simulations, these 
applications of the simulation results, as well as the mat- 
ter power spectrum and scalar field profile/evolution, will 
be presented in forthcoming works. 



Acknowledgments 

The authors thank Alexander Knebe, Claudio Llinares 
and Xufcn Wu for their helps in technical problems about 
the N-body code and its implementation, John Barrow, 
Carsten van der Bruck, Anne Davis, George Efstathiou, 
David F. Mota and Douglas Shaw for encouragements to 
finish this work and discussions in the process, and Luca 
Amendola for listening about the work and helpful in- 
formation at earlier stages. B. Li acknowledges supports 
from Overseas Research Studentship, Cambridge Over- 
seas Trust, DAMTP and Queens' College. We are also 
indebted to the HPC-Europa Transnational Access Visit 
programme for its support and the Lorentz Center and 
Leiden Observatory for hospitality when part of this work 
is undertaken. The N-body simulations are performed on 
the SARA supercomputer in the Netherlands. 



[1] E. J. Copeland, M. Sami and S. Tsujikawa, 

Int. J. Mod. Phys. D 15, 1753 (2006). 
[2] J. KhouryandA. Weltman, Phys. Rev. Lett. 93, 171104 

(2004). 

[3] J. Khoury and A. Weltman, Phys. Rev. D 69, 044026 
(2004). 

[4] D. F. Mota and D. J. Shaw, Phys. Rev. Lett. 97, 151102 
(2006). 



[5] D. F. Mota and D. J. Shaw, Phys. Rev. D 75, 063501 
(2007). 

[6] P. Brax, C. van de Bruck, A. -C. Davis, J. Khoury 
and A. Weltman, Phys. Rev. D 70, 123518 (2004). 

[7] S. M. Carroll, V. Duvvuri, M. Trodden and 
M. S Turner, Phys. Rev. D 70, 043528 (2004). 

[8] S. Nojiri and S. D. Odintsov, Phys. Rev. D 68, 123512 
(2003). 



17 



I. Navarro and K. van Acoleyen, JCAP 02, 022 (2007). 
B. Li and J. D. Barrow, Phys. Rev. D 75, 084010 (2007). 
W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007). 
P. Brax, C. van de Bruck, A. -C. Davis and D. J. Shaw, 
Phys. Rev. D 78, 104021 (2008). 

E. Bertschinger, Ann. Rev. Astron. Astrophys. 36, 599 

(1998) . 

E. V. Linder and A. Jenkins, Mon. Not. R. Astron. Soc. 
346, 573 (2003). 

R. Mainini, A. V. Maccio, S. A. Bonometto and 
A. Klypin, Astrophys. J. 599, 24 (2003). 
V. Springe! and G. R. Farrar, Mon. Not. R. Astron. Soc. 
380, 911 (2007). 

M. Kesden and M. Kamionkowski, Phys. Rev. Lett. 97, 
131303 (2006); Phys. Rev. D 74, 083007 (2006) 

G. R. Farrar and R. A. Rosen, Phys. Rev. Lett. 98, 
171302 (2007). 

J. A. Keselman, A. Nusser and P. J. E. Peebles (2007), 
|arXiv:0902.3452 [astro-ph]. 

A. V. Maccio, C. QuercelUni, R. Mainini, L. Amendola 
and S. A. Bonometto, Phys. Rev. D 69, 123516 (2004). 
M. Baldi, V. Pettorino, G. Robbers and V. Springe! 
(2008), afXiv108l2. 3901 [astro-ph]. 

P. Brax, C. van der Bruck, A. -C. Davis and 
A. M. Green, Phys. Lett. B 633, 441 (2006). 

H. Oyaizu, Phys. Rev. D 78, 123523 (2008). 

H. Oyaizu, M. Lima and W. Hu, Phys. Rev. D 78, 
123524 (2008). 

L. Amendoia, D. Polarski and S. Tsujikawa, 
Phys. Rev. Lett. 98, 131302 (2007). 

G. F. R. Ellis and M. Bruni, Phys. Rev. D 40, 1804 
(1989). 

A. Challinor and A. Lasenby, Astrophys. J. 513, 1 

(1999) . 

A. M. Lewis, A. Challinor and A. Lasenby, Astro- 
phys. J. 538, 473 (2000). 

S. Tsujikawa, Phys. Rev. D 76, 023514 (2007). 
A. Knebe, A. Green and J. Binney, Mon. Not. R. As- 
tron. Soc. 325, 845 (2001). 

A. Brandt, Math, of Comp. 31, 333 (1977). 

W. H. Press, S. A. Teukolsky, W. T. Vetterling and 

B. P. Flannery, Numerical Recipes in C. The Art of Sci- 
entific Computing (Cambridge University Press, Cam- 
bridge 1992), second ed. 

W. L. Briggs, V. E. Henson and S. F. McCormick, 

A Multigrid Tutorial (Society for Industrial and Applied 

Mathematics, Philadelphia 2000), second ed. 

Ya. B. Zel'dovich, Astron. Astrophys. 5, 84 (1970). 

G. Efstathiou, M. Davis, S. D. M. White and 

C. S. Frenk, Astrophys. J. Su ppl. 57, 241 (1985). 
E. Bertschinger 1970), arXiv: |astro^h/9506070 



C. Llinares, H. Zhao and A. Knebe, Astrophys. J. 695, 
L145 (2009). 



APPENDIX A: THE PERTURBATION 
EQUATIONS 



where Ua is the 4-velocity of an observer with respect to 
which 3-1-1 space-time splitting is made, hab = gab — UaUb 
the projection tensor used to obtain covariant tensors 
perpendicular to Ua- Here iTab is the projected symmet- 
ric tracefree anisotropic stress, q is the vector heat flux 
and p and p respectively the energy density and isotropic 
pressure. These quantities could be obtained from Tab 
through the relations 



P 

P 
qa 

T^ab 



TabU''u\ 



KhtT, 



cd 



Phab 



(A2) 



Based on these, the components of the stress energy ten- 
sor (up to first order) in this model is summarized in the 
following table: 



TABLE L The decomposition of energy momentum tensor in 
the 3 + 1 formalism. 



Matter 


P 


P 


qa 


T^ab 


7 


Pi 


zPi 


q^ia 




V 


Pv 




q^a 


T^uab 


Baryons 


PB 




qBa 





CDM 


PCDM 





qCDMa 





Coupled CDM 


C(</')/OCDM 





C{(p)qcDMa 







+ y(^) 


- F(^) 








Throughout this paper an overdot denotes the derivative 
with respect to the cosmic time t and is the covariant 
spatial derivative perpendicular to Ua (up to first order 
in perturbation). Note that it is the coupled CDM quan- 
tities which appear in the (background and perturbed) 
Einstein equations, and also note that we have included 
the massless neutrinos into the model here. 

The five constraint equations in the model are given as 





i^qa 

Bab 



T: H V CTah + V Wab, 



V'Bab 



2„ 



— K, 

2 



V TTab + 



^b)ec " 1 



2^ 



Vc^d + {p+p)-lUcd 



cd^b 



(A3) 
(A4) 
(A5) 
(A6) 
(A7) 



Consider the decomposition of the stress energy tensor 

Tab 



Here, eabcd is the covariant permutation tensor, £ab and 
Bab are respectively the electric and magnetic parts of 
the Weyl tensor Wabcd, given respectively through Sab 



Tab = TTab + 2q(^aUb) + pUaUb - phab (Al) u'^u'^Wacbd and Bab 



-u'^u^'^eac^yVefbd- 0,aab,rUab 
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come from the decomposition of the covariant derivative 
of 4-vclocity 

VaUb = (Jab + ^ab + + "a^b (AS) 

with A being the acceleration, 9 = V^Uc — 3d /a the 
expansion scalar, zuat — '^[a^b] and aab the shear. Note 
that 6 in the above section has a completely different 
meaning. 

In addition, the seven propagation equations are: 



Note that in all the above equations p, p, qa and Tiab are 
all the total quantities contributed by all matter species: 



p = P-I + Pu + PB + C{ip)pcT>M + 2"^^ + 

Qa = q'ya+qua+qBa+C{ifi)qcDMa + 'P^a'P, 
T^ab = T^-fab + T^uab- (A20) 



Qa 
Qa 



p+ip + p)9 + V''qa = 

Oqa + {p+p)Aa-WaP + ^''Trab = 0; 
Oqa + (P+ P)Aa ~ ^aP + V V^f, 



9± 
c 



(Jab + g^fafc 



V(aA) + ^ab 



1 
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T^ab H 


~ ^Ol^ab 












Lb + i 



■U3 + -Qw - V[aA] 



Bab + OBab + V'£d[aeb)ecU' + -KV%d(,eb)e/w' 



:(a;10) 

^(a;ii) 

:(a;12) 

■ im) 

:(M5) 



where the angle bracket means taking the trace-free part 
of a quantity. Note that the first of Eq. I|A10[) is for 
normal matter while the second is for the coupled dark 
matter. 

Besides the above equations, it is useful to express the 
projected Ricci scalar R into the hypersurfaces orthogo- 
nal to M° as 



R ^ 2kp- 



(A16) 



The spatial derivative of the projected Ricci scalar, r]a 
^aS/aR, is given as 



2a - 

T]a = anWaP - —OWat 



(A17) 



and its propagation equation 
29 2a 



9Vay Ab ~ aKVa^ qb- (A18) 



As we are considering a spatially flat universe, the spa- 
tial curvature must vanish on large scales which means 
that R = 0. Thus, from Eq. (jA16p . we can obtain the 
Friedman equation 



Kp. 



(A19) 



Finally, there is the perturbed scalar field EOM 



PCDM^ = 0(A21) 



APPENDIX B: EQUATIONS IN k SPACE 

To put the above perturbation equations into numeri- 
cal calculation, we need to write them in the fc-space. As 
we want to present the complete framework to study the 
structure formation in coupled scalar field models, here 
we also list those equations. 

First of all, the change into fc-space is accomplished 
with the aid of following Harmonic expansions: 

Xa=aVaP^Y.''^Qa qa ^ 

k k 

a 



^ab = J2 ^Qab 2a = a\/ a9 = ^ -ZQ^ 



k k 

ha = Vaa - V khQ^^ Aa = y] - AQ^ 

k k 

£-b = -Y. ^^Qab ^aV = E -^Qa (Bl) 

k k 



where — f'^aQ'' with Q'' being the zero order eigen- 
functions of the comoving Laplacian a^V^ {a^V^Q^ = 
fc^Q'') and = fV(a<5b)- ^o*^ ^^^^y consider 

scalar mode perturbations in the present paper and so 
shall neglect the second order quantities such as Bab and 

With these we list the equations that will be used in the 
numerical calculation of the linear large scale structure. 



19 



a. Background Evolution for Energy Densities 



Propagations of Heat Fluxes 



P7 + 4— 


= 0, 


(B2) 


a 


= 0, 


(B3) 


Pb + 3 — PB 


= 0, 


(B4) 


a' 

PCDM + 3 — PCDM 


= 0. 


(B5) 
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^ \ P^ 
3 V P7 



-aneCTT \v.y - -vb 



«B + — (1 - 3c2)ub + kA 
a 



-fee? A 



B + aUeCTT -VB 

PB \3 



v~ 



^CDM ^ ""CDM + kA+ —ip VCDM - 7T«^? 



= (BIO) 
= (B^1) 
= (B12) 

: (BI3) 



Note that because our definition of the dark matter en- 
ergy density includes only the mass density, but not the 
contribution from its coupling to the scalar field, so pcDM 
evolves as in ACDM. One can also understand this as fol- 
lows: the coupling to the scalar field produces a fifth force 
on the dark matter particles, making their trajectories be 
non-geodesic, but as we shall see below, the fifth force is 
spatial (perpendicular to the worldline of the dark mat- 
ter particle) and cannot change the mass of dark matter 
particles. Consequently the averaged dark matter mass 
density is the same as in ACDM. One can of course define 
the energy momentum tensor for dark matter as includ- 
ing C{(p), and in this case Eq. (jBSp will no longer be valid 
and we end up with varying mass dark matter particles. 



b. Propagations of Spatial Gradients of Densities 



d. Other Propagation Equations 



k[a' + —a] -k^(A + (j)) + -nlla^ = 0, (B14) 
^ a J 2 



a 
a 



1 

2 



n' H n - k{p + p)<7 - kq 

a 



= 0, (B15) 



kri' + 2—kA + Kqa^ = 0, (B16) 
a 



in which we have used 
P 



Pj + P^ + PB + C{ip)pcDM + ^f'^ + V{'flP17) 

^ ^ ^ ' (B18) 

(B19) 
(B20) 



P = -p, + -p. + —,v''-V{^l 

k 



n = n-, + n^ 



X = p^A^ + p^A^ + p^A^+ C{(p)PcoM'^CDM 

+ + f 4t^^" + ^^ + C^Pcdm) e (B21) 

and defined the density contrast Ai = Xi/ pi for matter 
species i. 



A' 



A' 



7 



4 a' 
A' + -kZ - 4—A + kv. 
' 3 a 

4 a' 
AI + -kZ - 4:—A + kv^ 
6 a 

kZ - 3— A + kvS] + 3— c^Ab 
a la 



e. Constraint Equations 



CDM 



kZ — i—A + kvcDM 
a 



0, (B6) 

0, (B7) 

0, (B8) 

0. (B9) 



k-^fj} + -K 



k'^{Z - fj) + -Kqa'^ 



k(n + X) + 3—q 
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k^ri + 2—kZ - nXa^ 
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(B22) 


0, 


(B23) 


0. 


(B24) 
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f. Scalar Field Equations of Motion 



a' , dV{v) 2 2 

2— H a +PCDM — a 

a O'f dip 



+ 2ip" + —(p' A + a^C^pcDMAcDM 



Instead, we can define a new variable u according to 

e" + l = eX. (CI) 

By this, u still takes value in (—00,00), e" G (0, 00) and 
= (1P25) thus e (1, 00) which ensures that x is positive definite 
in numerical solutions. Besides, e^^ = [1 + e"]'^ so that 
there will be no exponential-of-exponential terms, and 
the only exponential is what we have for the potential 
itself. P = —1 above. 

Then the Poisson equation becomes 



+ {kZ + A')ip' ^ (B26) 



g. The Friedmann Equations 



a 
a 



31^ 

a 



npa 



6 



(B27) 
(B28) 



When it comes to perturbation calculations, we need 
to fix a gauge, i.e., choose a Ua- One possibility is to use 
the 4- velocity of dark matter particles as our Ua- In this 
case the dark matter heat flux is zero and according to 
Eq. (jB13[) we will simply have A — and this relation 
can then be used to replace the A^s appearing in all the 
above equations. Another possibility is to choose Ua such 
that A — 0: in this case vcdm will become nonzero and 
we need to dynamically evolve it. Below in the numerical 
calculations we shall adopt the second possibility. 



APPENDIX C: DISCRETIZED EQUATIONS FOR 
N-BODY SIMULATIONS NON-LINEAR REGIME 



In the MLAPM code the partial differential equation 
Eq. (pS)) is (and in our modified code Eq. P5)) will also 
be) solved on discretized grid points, and as such we must 
develop the discretized versions of Eqs. ([^ - to be 
implemented into the code. 

But before going on to the discretization, we need to 
address a technical issue. As the potential is highly non- 
linear, in the high density regime the value of the scalar 
field ^/K^p will be very close to 0, and this is potentially 
a disaster as during the numerical solution process the 
value of \/Kip might easily go into the forbidden region 
(f < [23;]. One way of solving this problem is to define 
Y = xe" in which x is the background value of x, as in 
[23|. Then the new variable u takes value in (—00, 00) so 
that e" is positive definite which ensures that x > 0. 
However, since there are already exponentials of x ij^i 
the potential, this substitution will result terms involving 
exp [exp(u)], which could potentially magnify any numer- 
ical error in u. 



-n 



CDM 



1- (1-t-e")'^ 



(C2) 



where we have defined fly = i^V{'p)/'^H^ which is deter- 
mined by background cosmology, the quantity e'^'t' is 
also determined solely by background cosmology. These 
background quantities should not bother us here. 
The scalar field EOM becomes 



rV- 



-Vu 



37^^CDMPc (1 + e^y + 



3Ai/3f)voaMl + e")'^ 



1 - (1 + e"^ f 



-37f2cDMe 



3/i/3r2voa3e'3^^ 



1^1 _ g/3\/K¥'j 



(C3) 



in which we have used the fact that x = log(l + e") 

= i-^e" ^^'^ moved all terms depending only on 
background cosmology (the source terms) to the right 
hand side. 

So, in terms of the new variable w, the set of equations 
used in the N-body code should be 



dpc 



2 

c^7 



1 



-Vu 



(C4) 
(C5) 



plus Eqs. (102) IC3|) . These equations will ultimately be 
used in the code. Among them, Eqs. (jC2|IC5P will use the 
value of u while Eq. (|C3p solves for u. In order that these 
equations can be integrated into MLAPM, we need to 
discretize Eq. (jC3|) for the application of Newton-Gauss- 
Seidel iterations. 

To discretize Eq. (IC3|, let us define b = The 
discretization involves writing down a discretion version 
of this equation on a uniform grid with grid spacing h. 
Suppose we require second order precision as is in the 
standard Poisson solver of MLAPM, then Vu in one di- 
mension can be written as 



2h 



(C6) 
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where a subscript j means that the quantity is evaluated The factor h in V- (6Vu) makes this a standard variable 

on the j-th point. Of course the generalization to three coefficient problem. We need also discretize 6, and do it 

dimensions is straightforward. in this way (again for one dimension): 
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where we have defined bj_^_i = {bj + bj+i) /2 and bj_i = 
{bj-i + bj) /2. This can be easily generalize to three di- 



mensions as 



J 



(C7) 



V • (bVu) 



1 








1 






1 







Then the discrete version of Eq. (|C3P is 
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Then the Newton-Gauss-Seidel iteration says that we can rate) solution uf^i, as 
obtain a new (and often more accurate) solution of u, 
u^Yk^ using our knowledge about the old (and less accu- 



i,j,k 



old 

ijj^k 
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The old solution will be replaced by the new solution to Gauss-Seidel sweeping scheme. Note that 
once the new solution is ready, using the red-black 
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In principle, if we start from a high redshift, then the 
initial guess of Uij^k could be such that the initial value 
of X ill 8-11 space is equal to the background value Xi 
because anyway at this time we expect this to be approx- 
imately true. For subsequent time steps we could use the 
solution for Uij^k at the previous time step as our initial 
guess; if the time step is small enough then we do not ex- 



pect u to change significantly between consecutive times 
so that such a guess will be good enough for the iteration 
to converge fast. 

In practice, however, due to specific features and al- 
gorithm of the MLAPM code [30| , the above procedure 
may be slightly different in details. 



